function vamdl = fitvadist(vadata,numcomp,emopt)

if nargin<3, emopt = {}; end

%-- Expand Data --%;
[numTCH,num_outcome] = size(vadata);

if num_outcome==1
    Y = {vadata.Y};
    Z = {vadata.Z};
    ZF = {vadata.ZF};
else
    assert(isequal(range(reshape([vadata.teacherid],size(vadata)),2),zeros(numTCH,1)))
    Y = arrayfun(@(j) cat(1,vadata(j,:).Y),1:numTCH,'unif',0);
    Z = arrayfun(@(j) blkdiag(vadata(j,:).Z),1:numTCH,'unif',0);
    ZF = arrayfun(@(j) blkdiag(vadata(j,:).ZF),1:numTCH,'unif',0);
end

if ~(rcond(cat(1,Z{:})'*cat(1,Z{:}))>eps) 
    error('match variables are not full rank') 
end

numvars = size(cat(1,Z{:}),2);

%-- Generate fixed effect data --%;
disp('Generate fixed effect data')

tchfe = cellfun(@(zf,y) (zf'*zf) \ (zf'*y),ZF,Y,'unif',0);
tchfeMAP = cellfun(@(zf,z) (zf'*zf) \ (zf'*z),ZF,Z,'unif',0);

if num_outcome==1
    resid = cellfun(@(y,zf,fe) y - zf*fe,Y,ZF,tchfe,'unif',0);
    resid = cat(1,resid{:});
    df = numel(resid) - numel(cat(1,tchfe{:}));
    residvar = sum(resid.^2)/df;
    tchfeERRCOV = cellfun(@(zf) residvar*((zf'*zf) \ eye(size(zf,2))),ZF,'unif',0);
    adjR2 = 1 - (sum(resid.^2)/df)/var(cat(1,Y{:}),0,1);
else
    if ~all(cellfun(@(i) all(diff(i)>0),{vadata.sid}))
        error('sid data is out of order')
    end
    I = cell([1 numTCH]);
    resid = cell([1 numTCH]);
    for j = 1:numTCH
        uni_sid = unique(cat(1,vadata(j,:).sid));
        sid_ind = cellfun(@(sid) ismember(uni_sid,sid),{vadata(j,:).sid},'unif',0);
        I{j} = cat(2,sid_ind{:});
        resid{j} = zeros(size(I{j}));    
        resid{j}(I{j}) = Y{j} - ZF{j}*tchfe{j};
    end
    resid = cat(1,resid{:});
    df = sum(cat(1,I{:})) - sum(reshape(cellfun(@(x) size(x,2),{vadata.ZF}),size(vadata)));
    df = diag(df);
    for j = 1:numTCH
        MZ = cellfun(@(z) eye(size(z,1)) - z*((z'*z) \ z'),{vadata(j,:).ZF},'unif',0);
        for co = 1:num_outcome
            for ro = (co+1):num_outcome
                OV = I{j}(:,ro) & I{j}(:,co);
                A = eye(size(I{j},1));
                A = A(I{j}(:,ro),I{j}(:,co));
                df(ro,co) = df(ro,co) + trace(MZ{ro}(OV(I{j}(:,ro)),:)'*MZ{co}(OV(I{j}(:,co)),:)*A');
            end
        end
    end
    residcov = (resid'*resid)./df;
    residcov(triu(true(num_outcome),1)) = residcov(tril(true(num_outcome),-1))';
    tchfeERRCOV = cell([1 numTCH]);
    for j = 1:numTCH
        ERRCOV = kron(residcov,speye(size(I{j},1)));
        ERRCOV = ERRCOV(I{j},I{j});
        zzfINV = (ZF{j}'*ZF{j}) \ eye(size(ZF{j},2));
        tchfeERRCOV{j} = zzfINV*ZF{j}'*ERRCOV*ZF{j}*zzfINV;
    end
    adjR2 = nan;
end
tchfeERRCOV = cellfun(@(x) (x + x')./2,tchfeERRCOV,'unif',0);

%-- Starting Values --%;
disp('Getting Starting Values')
ucz2 = mean(cat(1,Z{:}).^2);
ls_hat = cellfun(@(z,y) (z'*z + numvars*diag(ucz2)) \ z'*y,Z,Y,'unif',0);
ls_hat = cat(2,ls_hat{:})';
ls_hat(ls_hat==0) = nan;
ls_E = mean(ls_hat,'omitnan')';
ls_C = cov(ls_hat,"partialrows");
vahat = cellfun(@(Y,H,Q) mvnfilter(Y,H,ls_E,ls_C,Q),tchfe,tchfeMAP,tchfeERRCOV,'unif',0);
vahat = cat(2,vahat{:})';

vadist = struct('mean',zeros([numvars numcomp])...
                ,'cov',zeros([numvars numvars numcomp])...
                ,'pr',zeros([1 numcomp]));
idx = kmeans(vahat,numcomp,'Replicates',100,'MaxIter',1000);
for c = 1:numcomp
    vadist.mean(:,c) = mean(vahat(idx==c,:));
    vadist.cov(:,:,c) = diag(var(ls_hat(idx==c,:),0,'omitnan'));
    vadist.pr(c) = mean(idx==c);
end

%-- Estimation --%;
%chkderiv(tchfe,tchfeMAP,tchfeERRCOV,vadist);
[vadist,output] = emqnalgo(vadist,tchfe,tchfeMAP,tchfeERRCOV,emopt{:});
if num_outcome==1
    vadist.residvar = residvar;
else
    vadist.residcov = residcov;
end

%-- Standard Errors --%;
[~,~,~,SCR] = emstep(vadist,tchfe,tchfeMAP,tchfeERRCOV);
VCOV = (SCR'*SCR) \ eye(size(SCR,2));

p2p = @(x) vecparm(x,numvars,numcomp);
parms = p2p(vadist);
calcSE = @(fun) sqrt(deltamethod(fun,parms,VCOV));

vadistSE = p2p(sqrt(diag(VCOV)));
vadistSE.pr = calcSE(@(x) p2p(x).pr);

vamom = gmmom(vadist);
vamomSE = cellfun(@(f) calcSE(@(x) gmmom(p2p(x)).(f)),fieldnames(vamom),'unif',0);
vamomSE = cell2struct(vamomSE,fieldnames(vamom));

%-- Organize Output --%;
vamdl = struct('VAdistEstimates',vadist ...
    ,'VAdistStandardErrors',vadistSE ...
    ,'VAmomentEstimates',vamom ...
    ,'VAmomentStandardErrors',vamomSE ...
    ,'NumParameters',output.numparms ...
    ,'NumRandomCoefficients',numvars ...
    ,'NumMixtureComponents',numcomp ...
    ,'NumObservationalUnits',numTCH ...
    ,'OptimizationInfo',rmfield(output,'numparms'));

vaout = [num2cell([vadata(:,1).teacherid]') cell([numTCH 3])];
for j = 1:numTCH
    [vaout{j,2:end}] = vaposterior(vadist,vadata(j,:));
end
vaout = cell2struct(vaout',{'teacherid','vahat','vadist','ll'});

ll = sum([vaout.ll]);    
vamdl.LogLikelihood = ll;
vamdl.AIC = 2*vamdl.NumParameters - 2*ll;
vamdl.BIC = log(numTCH)*vamdl.NumParameters - 2*ll;
vamdl.RsquaredAdjusted = adjR2;

vamdl.ValueAddedEstimates = vaout;

[~,fe_idx] = cellfun(@(W) ismember(round(W,10),eye(numvars),'rows'),tchfeMAP,'unif',0);
lsout = struct('teacherid',num2cell([vadata(:,1).teacherid]) ...
    ,'coefid',cellfun(@(idx) nonzeros(idx),fe_idx,'unif',0) ...
    ,'coef',cellfun(@(fe,idx) fe(idx~=0,:),tchfe,fe_idx,'unif',0) ...
    ,'coefcov',cellfun(@(V,idx) V(idx~=0,idx~=0),tchfeERRCOV,fe_idx,'unif',0))';

vamdl.LeastSquaresEstimates = lsout;

end
